1 Load Data

OD_data <- read_excel(file.path("data", "2018_Silverman_et_al-OD_data.xlsx"))

2 Growth curves

2.1 Bootstrapped means for OD data

n_bootstrap <- 1000
OD_summary <- OD_data %>% 
  nest(-organism, -pN2, -type, -day) %>% 
  mutate(
    # number of samples
    n_samples = map_int(data, nrow),
    # bootstrap data to estimate mean and standard error
    bootstrap = map(data, function(.x) {
      if (nrow(.x) > 1) boot(data = .x$OD750, statistic = function(x, idx) mean(x[idx]), R = n_bootstrap)
      else list(t = .x$OD750)
    }),
    # bootstrapped based mean and standard error
    OD750_mean = map_dbl(bootstrap, ~mean(.x$t)),
    OD750_mean_se = map_dbl(bootstrap, ~sd(.x$t)),
    # error bars upper and lower
    OD750_mean_ue = OD750_mean + OD750_mean_se,
    OD750_mean_le = OD750_mean - OD750_mean_se,
    # pN2 formatted for legend
    pN2_formatted = sprintf("%0.1f bar", pN2)
  )

2.2 Visualization of all data and means

# plot OD both on linear and log scale
bind_rows(
  # linear scale data
  mutate(OD_summary, scale = "OD750"),
  # log scale data
  mutate(filter(OD_summary, OD750_mean > 0, day > 0, type == "sample"), scale = "log10 (OD750)", 
         data = map(data, ~mutate(.x, OD750 = log10(OD750))),
         OD750_mean = log10(OD750_mean),
         OD750_mean_le = log10(OD750_mean_le),
         OD750_mean_ue = log10(OD750_mean_ue))
) %>% 
  ggplot() +
  # plot-wide aesthetics
  aes(x = day, y = OD750_mean, fill = pN2_formatted, color = pN2_formatted, shape = type) +
  # all data (small translucent points)
  geom_point(data = function(df) unnest(df, data), map = aes(y = OD750), size = 2, alpha = 0.5) +
  # error bars (1 SE) for averages
  geom_errorbar(data = function(df) filter(df, !is.na(OD750_mean_se)), 
                map = aes(ymin = OD750_mean_le, ymax = OD750_mean_ue), width = 0.2, alpha = 0.6) +
  # averages and lines connection them
  geom_line(size = 1) +
  geom_point(size = 4, color = "black") +
  # OD cutoff indicator for growth rate calculations
  geom_hline(
    data = data_frame(scale = c("OD750", "log10 (OD750)"), cutoff = c(0.4, log10(0.4))),
    map = aes(yintercept = cutoff, color = NULL, fill = NULL, shape = NULL)
  ) +
  # panels
  facet_grid(scale~organism, space = "free_x", scales = "free") + 
  # scale
  scale_x_continuous(breaks = 0:25) +
  scale_color_manual(values = cbPalette) +
  scale_fill_manual(values = cbPalette) +
  scale_shape_manual(values = c(21:26)) +
  # labels
  labs(x = "day", y = TeX("$OD_{750}$"), color = expression("pN"[2]), fill = expression("pN"[2]), size = 5) +
  # theme
  theme_figure(grid = FALSE, text_size = 20) +
  guides(fill = guide_legend(override.aes = list(shape = 21)))

3 Growth rates

3.1 Analysis

# calculate growth rates
growth_rates <- 
  OD_data %>% 
  filter(day > 0, OD750 > 0, type == "sample", pN2 > 0, OD750 > 0.01, OD750 <= 0.4) %>% 
  nest(-organism, -pN2) %>%
  mutate(
    # bootstrap slope from regression model
    bootstrap = map(data, ~boot(data = .x, R = n_bootstrap,
      statistic = function (x, idx) {
        # model fit for resampled data
        m <- lm(log(OD750/min(OD750)) ~ day, data = x[idx,])
        r2 <- glance(m)$r.squared
        slope <- tidy(m)$estimate[2]
        return(c(slope, r2))
      })),
    growth_rate.day = map_dbl(bootstrap, ~mean(.x$t[,1])),
    growth_rate_se.day = map_dbl(bootstrap, ~sd(.x$t[,1])),
    r2 = map_dbl(bootstrap, ~mean(.x$t[,2]))
  ) %>% 
  select(-bootstrap, -data)
growth_rates %>% write.xlsx(file.path("data", "2018_Silverman_et_al-growth_rate_data.xlsx"))
growth_rates

3.2 Visualization

# plot growth rates
growth_rates %>% 
  ggplot() + 
  aes(pN2, growth_rate.day, color = organism, shape = organism) +
  geom_errorbar(aes(ymin = growth_rate.day - growth_rate_se.day, ymax = growth_rate.day + growth_rate_se.day), width = 0) +
  geom_point(size = 5) +
  scale_color_manual("Species", values = cbPalette) +
  scale_shape_manual("Species", values = c(16, 15)) + 
  labs(y = "growth rate [1/day]", x = TeX("$pN_{2}$ \\[bar\\]")) +
  theme_figure()